Contents

Surface Matching: example on human scapula

1.1. Surface Matching: example on human scapula

This notebook presents a method to match a point cloud taken with any scanning system to a digital 3D scan in stl format.

1.1.1. Introduction

1.1.1.1. Method outlines

Surface matching is based on an optimization strategy. The objective of this optimization is to identify the parameters of the reference change allowing to pass from the coordinate system of the STL file to the coordinate system of the scan system.

This change of reference frame is formalized by a passage matrix formed by a rotation vector and a translation vector according to the Rodrigues formalism. There are thus a total of 6 DOF.

The optimization parameters of the cost function will therefore be these two vectors. The residual is defined by the signed distance between the points coming from the measurement systems after the change of reference frame and the stl mesh. This distance is calculated between the points of the closest scanning system and the STL mesh.

This problem being badly formulated, it is necessary to establish an initial solution in order to have a correct result.

This is why in a first step, we establish an initial using points that are identified on the STL model but also on the real model. A minimum of 3 points is necessary to block the 6DOF.

Then, in a second step, an optimization will be performed with the whole scan taking as initial solution the one previously calculated.

1.1.1.2. Modules importations

import trimesh
import numpy as np
from stl import mesh
from scipy import optimize
import plotly.graph_objects as go
import pandas as pd
import gbu

1.1.1.3. Case of study

Left Human scapula stl file

def stl2mesh3d(stl_mesh):
    # stl_mesh is read by nympy-stl from a stl file; it is  an array of faces/triangles (i.e. three 3d points)
    # this function extracts the unique vertices and the lists I, J, K to define a Plotly mesh3d
    p, q, r = stl_mesh.vectors.shape  # (p, 3, 3)
    # the array stl_mesh.vectors.reshape(p*q, r) can contain multiple copies of the same vertex;
    # extract unique vertices from all mesh triangles
    vertices, ixr = np.unique(
        stl_mesh.vectors.reshape(p * q, r), return_inverse=True, axis=0
    )
    I = np.take(ixr, [3 * k for k in range(p)])
    J = np.take(ixr, [3 * k + 1 for k in range(p)])
    K = np.take(ixr, [3 * k + 2 for k in range(p)])
    return vertices, I, J, K


def create_mesh3D(
    stl_file="scapula.stl",
    title="Mesh3d from a STL file<br>AT&T building",
    colorscale=None,
):

    if colorscale is None:
        colorscale = [[0, "#e5dee5"], [1, "#e5dee5"]]
    vertices, I, J, K = stl2mesh3d(mesh.Mesh.from_file(stl_file))
    x, y, z = vertices.T
    mesh3D = go.Mesh3d(
        x=x,
        y=y,
        z=z,
        i=I,
        j=J,
        k=K,
        flatshading=True,
        colorscale=colorscale,
        intensity=z,
        name="AT&T",
        showscale=False,
    )

    layout = go.Layout(
        paper_bgcolor="rgb(1,1,1)",
        title_text=title,
        title_x=0.5,
        font_color="white",
        width=800,
        height=800,
        scene_camera=dict(eye=dict(x=1.25, y=-1.25, z=1)),
        scene_xaxis_visible=False,
        scene_yaxis_visible=False,
        scene_zaxis_visible=False,
    )
    fig = go.Figure(data=[mesh3D], layout=layout)

    fig.data[0].update(
        lighting=dict(
            ambient=0.18,
            diffuse=1,
            fresnel=0.1,
            specular=1,
            roughness=0.1,
            facenormalsepsilon=0,
        )
    )
    fig.data[0].update(lightposition=dict(x=3000, y=3000, z=10000))

    return fig


def add_points(points, fig, name, scale=1e3, *args, **kwargs):
    fig.add_trace(
        go.Scatter3d(
            x=np.concatenate([points, points]).reshape(-1, 3)[:, 0] * scale,
            y=np.concatenate([points, points]).reshape(-1, 3)[:, 1] * scale,
            z=np.concatenate([points, points]).reshape(-1, 3)[:, 2] * scale,
            name=name,
            mode="markers",
            **kwargs
        )
    )
    return fig


fig = create_mesh3D(stl_file="scapula.stl", title="Scapula Left <br> stl file")
fig.show()

1.1.2. STL file

We define reference points on the STl that can be easily found on the real model.
A point on the glenoid, coracoid and acromion is selected.

p_glene_stl = np.array([0.02515289854664963, 0.03122873596800735, 0.036798809060995745])
p_coracoide_stl = np.array(
    [0.026787610816186087, 0.05609095153197788, -0.0043762069034760384]
)
p_acromion_stl = np.array(
    [-0.008420164259823573, 0.07140986464971243, 0.018059532550748165]
)
fig = add_points(points=p_glene_stl, fig=fig, name="glène")
fig = add_points(points=p_coracoide_stl, fig=fig, name="coracoide")
fig = add_points(points=p_acromion_stl, fig=fig, name="acromion")
fig.show()

The same points are easily found on the real model.
Here, the points are marked in green.

../_images/scapula.jpg

1.1.3. System scan device

1.1.3.1. Get scanned points from our system scan device

We get the points recorded by our system system scan device. Each reference points are labeled respectively and so the full scan of glenoid.

df_scan = pd.read_pickle("data_scan.p")
df_scan
p_glene_scapula p_acromion_scapula p_coracoide_scapula p_scan_scapula
p3d p3d p3d p3d
x y z x y z x y z x y z
0 0.014046 0.010801 0.118422 -0.023197 -0.022039 0.086035 -0.002085 -0.037214 0.123242 0.015003 -0.004692 0.11943
1 0.014152 0.01073 0.118276 -0.023097 -0.022086 0.085978 -0.002091 -0.037263 0.123189 0.015018 -0.004673 0.119397
2 0.014185 0.010854 0.11838 -0.023238 -0.022144 0.086024 -0.002113 -0.037188 0.123366 0.015072 -0.004562 0.119343
3 0.014362 0.01094 0.118554 -0.023217 -0.022149 0.086317 -0.001946 -0.037163 0.123109 0.015048 -0.004463 0.119347
4 0.014379 0.010899 0.118494 -0.023224 -0.022023 0.086214 -0.001959 -0.036706 0.12319 0.015013 -0.004507 0.11929
... ... ... ... ... ... ... ... ... ... ... ... ...
79 NaN NaN NaN NaN NaN NaN NaN NaN NaN 0.015566 -0.001602 0.126646
80 NaN NaN NaN NaN NaN NaN NaN NaN NaN 0.01562 -0.001514 0.126536
81 NaN NaN NaN NaN NaN NaN NaN NaN NaN 0.015541 -0.001624 0.126567
82 NaN NaN NaN NaN NaN NaN NaN NaN NaN 0.015546 -0.001497 0.126532
83 NaN NaN NaN NaN NaN NaN NaN NaN NaN 0.015578 -0.001497 0.126524

84 rows × 12 columns

1.1.3.2. Plot scanned points in their reference frame

Here the points are expressed in scapula reference frame

cols = set([cols[0] for cols in df_scan.columns])
fig = go.Figure()
for i, c in enumerate(cols):
    fig.add_trace(
        go.Scatter3d(
            x=df_scan[c].p3d.x,
            y=df_scan[c].p3d.y,
            z=df_scan[c].p3d.z,
            name=c,
            marker_size=2,
            mode="markers",
        )
    )
fig.show()

1.1.4. Create initial guess

1.1.4.1. Gathering input data of reference points

p_glene_scan = np.array(df_scan["p_glene_scapula"].p3d.mean())
p_coracoide_scan = np.array(df_scan["p_coracoide_scapula"].p3d.mean())
p_acromion_scan = np.array(df_scan["p_acromion_scapula"].p3d.mean())

tri_point_stl = np.concatenate([p_glene_stl, p_coracoide_stl, p_acromion_stl]).reshape(
    -1, 3
)
tri_point_scan = np.concatenate(
    [p_glene_scan, p_coracoide_scan, p_acromion_scan]
).reshape(-1, 3)

1.1.4.2. Compute initial guess

def unpack(X):
    return X.reshape(2, 3)


def cost(X):
    Rsc2stl, Tsc2stl = unpack(X)
    tri_point_out = gbu.utils.apply_RT(tri_point_scan, Rsc2stl, Tsc2stl)
    dist = abs(tri_point_stl - tri_point_out)
    pen_pin_radius = 0
    res = dist - pen_pin_radius

    return res.flatten()


X0 = np.zeros(6)
sol = optimize.least_squares(
    cost,
    X0,
    method="lm",
    ftol=1.0e-12,
    xtol=1.0e-12,
    gtol=1.0e-10,
)
X_pre = sol.x
R_inital_guess, T_initial_guess = unpack(X_pre)
print(f"Initial guess solution: rvec={R_inital_guess}, tvec={T_initial_guess}")

Initial guess solution: rvec=[ 2.23440972 -0.14576349  0.55625112], tvec=[-0.02428643  0.12368749  0.09225085]

1.1.5. Surface match with full scan

1.1.5.1. Gathering input data of full scan

p_full_scan_sc = np.array(df_scan["p_scan_scapula"])

# load full scapula expand
scale = 1e-3  # convert to meter
mesh_scapula = trimesh.exchange.load.load("scapula.stl", type="stl")
mesh_scapula.vertices *= scale

1.1.5.2. Optimization

def cost(X, p3d, mesh):
    Rsc2stl, Tsc2stl = unpack(X)
    p_stl = gbu.utils.apply_RT(p3d, Rsc2stl, Tsc2stl)
    dist = trimesh.proximity.signed_distance(mesh, p_stl)
    pen_pin_radius = 0
    res = dist - pen_pin_radius
    return res


sol = optimize.least_squares(
    cost,
    X_pre,
    method="lm",
    ftol=1.0e-12,
    xtol=1.0e-12,
    gtol=1.0e-10,
    args=(p_full_scan_sc, mesh_scapula),
)

R_sc2stl_sol, T_sc2stl_sol = unpack(sol.x)
print(f"Final solution: rvec={R_sc2stl_sol}, tvec={T_sc2stl_sol}")

Final solution: rvec=[ 2.2283109  -0.16742782  0.55011184], tvec=[-0.02341616  0.12482507  0.09038642]

1.1.5.3. Plot results

p_full_scan_stl = gbu.utils.apply_RT(p_full_scan_sc, R_sc2stl_sol, T_sc2stl_sol)

fig = create_mesh3D(stl_file="scapula.stl", title="Scapula Left <br> stl file")

fig = add_points(points=p_full_scan_stl, fig=fig, name="scanned_points", marker_size=10)

fig.show()

previous

1. Surface Matching